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Abstract. - We present a non- neutral stochastic model for the dynamics taking place in a meta- 
community ecosystems in presence of migration. The model provides a framework for describing 
the emergence of multiple ecological scenarios and behaves in two extreme limits either as the 
unified neutral theory of biodiversity or as the Bak-Sneppen model. Interestingly, the model 
shows a condensation phase transition where one species becomes the dominant one, the diversity 
in the ecosystems is strongly reduced and the ecosystem is non-stationary. This phase transition 
extend the principle of competitive exclusion to open ecosystems and might be relevant for the 
study of the impact of invasive species in native ecologies. 



Understanding the basic principles responsible for the 
biodiversity in ecosystems is a fundamental questions of 
scientific and social interest. In fact many ecosystems re- 
cently, undergo a loss of biodiversity [1,2] due to human 
activities or to non-native invasive species that in some 
cases might become dominant in the ecosystem. 

Recently the unified neutral theory of biodiversity and 
biogeography [3] has been shown to reproduce species abun- 
dance distributions [4, 5] of meta-communities and local 
communities and the distribution of lifetimes [6] that are 
found in ecological systems. The proposed theoretical 
framework is based on the assumption that all species 
are equivalent and assumes that species abundance dis- 
tribution is affected mainly by the drift resulting from a 
stochastic process of births, deaths and speciation. This 
theory raised some controversies in the field [7,8], and 
there is an increasing need of a theory which includes 
the stochastic features of the neutral model but also ac- 
counts for the competitive advantage that some species 
might have in their niche. Moreover a non-neutral theory 
of meta-communities biodiversity has the potential to in- 
clude also complex prey-predator foodweb relations [9-11] 
for which an increasing number of data is now available. 

In the physics community non-neutral ecological mod- 
els have been widely studied [11-13]. In this context a 
reference model is the Bak-Sneppen model [14,15], which 
describes the extinction process in a non-neutral ecosys- 
tem. In the model, a negative selection process takes place 



and at each step the less fit species of the ecosystem is re- 
moved, and mutations affect the species dependent on the 
removed one. This model leads to punctuated equilibrium 
and power-law extinction patterns found in fossil records. 

In this paper we show that these two paradigms (the 
neutral model and the Bak-Sneppen model) can be re- 
covered by a simple model of metacommunities ecological 
systems [11, 12] evolving following the principles of evo- 
lutionary dynamics, which include the fundamental in- 
gredients of reproduction, negative selection and migra- 
tion/speciation. In our stochastic model we assume that 
species are not equivalent but they are characterized by 
a quenched parameter that encodes for the adaptation of 
the species to their niche. We call this parameter the fit- 
ness of the species. High fitness will imply a competitive 
advantage of the species as in evolutionary models. In 
particular high fitness will increases the probability of re- 
production of a species, and reduce the death probability 
due to a negative selection process. In the model spe- 
ciations processes are included and we assume that they 
occur mainly as a consequence of migration events, in or- 
der to describe invasive species arriving at fixed rate in 
an ecosystem. This model reduces to the neutral theory 
of biodiversity in the limit in which all species have the 
same fitness. On the other side when we include the effect 
of evolutionary selection on species of heterogeneous fit- 
ness the model has a rich phenomenology. In the limit in 
which only negative selection is present, and considering 
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an extremal dynamics, we recover a dynamics that can be 
mapped into the mean-field Bak-Sncppen model [14,15]. 
Moreover the model shows a phase transition to a non- 
stationary state in which one species becomes the domi- 
nant species in the ecosystem and the biodiversity of the 
ecosystem is strongly reduced. This phase transition is a 
condensation phase transition that in some limiting cases 
can be mapped to a Bose-Einstein condensation and is 
closely related to the condensation transitions found in 
different evolving models [16-20]. Below this phase tran- 
sition, relevant dynamical effects take place and the con- 
densed species with a finite population change in time in 
a nontrivial way. 

Non-neutral birth-death dynamics of the ecosystem - We 
generalize the stochastic birth-death dynamics giving rise 
to the distribution of species in meta-communities in the 
neutral theory of biodiversity to a birth-death dynamics 
depending on quenched variables assigned to the species 
and that correspond to their niche adaptation. We call 
this quenched variable the fitness / of a species. We con- 
sider an ecosystem formed by N individuals belonging to 
S species with the number of individuals N fixed but a 
variable number of species S. Each individual belongs 
to a species i — 1, 2, . . . , S and each species i has fitness 
fi > 0. The birth-death dynamics is inspired by evolu- 
tionary dynamics and three different birth-death processes 
take place. We call this three processes replication, neg- 
ative selection and migration/speciation. Starting from 
a random initial condition we define the following birth- 
death process: 

• Replication - At rate a an individual of the species i 
with fitness fi is chosen to replicate with probability 
w(fi), where w(f) is an increasing function of the 
fitness, and is substituted by a random individual of 
the ecosystem. Therefore the resulting process is 



rii + 1 and nj — > nj — 1 



(1) 



and the species i and j are selected according to the 
probabilities Ilf and Ilj* respectively, defined as 



nf = and Ilf = ^. (2) 



This is the usual process of reproduction in popu- 
lation dynamics in which w(fi) describes the mean 
number of expected off-springs of the individual i. 

• Negative selection - At rate 7 an individual of a 
species i with fitness fi is removed from the popu- 
lation with probability B(fi), which is a decreasing 
function of the fitness fi and is substituted by an off- 
spring of a random individual. The resulting process 
is 



and the species i and j are selected according to the 
probabilities U^ s and nj^ 5 respectively, defined as 
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and 



X (4) 



This is a term of negative selection and describes the 
struggle for survival of the individuals. In particular 
B(fi) we take 



B(fi) = l/w(fi) 



(5) 



describing the Darwin concept of removal of the less 
fit species. More general form of B(f) might instead 
represent barrier between fitness peaks in some com- 
plex fitness landscape. 

This type of process is the pivotal process of the Bak- 
Sneppen model which is in addition also extremal (at 
each time-step the less fit species is removed from the 
ecosystem or migrates away from the metacommunity 
). We observe here that in our model when a = 0, 
7 < 1 and (3 — ► 00, we recover an extremal dynamics 
that is reminiscent of the Bak-Sneppen model [14, 15]. 

One might naively think that the effect of the negative 
selection process and of the reproduction process is 
equivalent. This is not the case in general as it will 
be shown in the following. 

• Migration/Speciation At rate /j, = 1 a migra- 
tion/speciation occurs and a random individual is 
substituted with an individual of a new species i = 
S + 1. The new species is assigned a fitness fs+i 
randomly drawn from a fixed distribution p(f). The 
resulting process is 



ns+i — > 1 and 



nj - 1 



(6) 



M 



where the species j is selected with probability II 

and the fitness fs+i is drawn from a distribu- 
tion /?(/). We consider here on purpose only the case 
in which the new species have fitnesses independent 
on the fitnesses of the species present in the ecosys- 
tem. This is the case of speciations due to migration 
from another metacommunity which is relevant for 
the study of invasive species in given ecosystems. 

In particular we consider the case in which the function 

w(fi) = e^. (7) 

where [3 is a tunable parameter. The assumption (7) is by 
no means a limitation since the model could be defined 
equivalently only in terms of the distribution of the vari- 
ables {w}'s. Moreover the choice of the assumption (7) 
has the following advantages: 



1 and n,- 



1 



(3) 



i) It ensure that we have always w(fi) > 0. 
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• ii) It allows to mimic the impact of the environment 
in driving the system out of neutrality, by changing 
(3. In fact (3 tunes the relevance of the fitness in the 
stochastic dynamics. In particular when (3 — we 
have Wi = 1 Vi and in the dynamical process all 
species are equivalent, recovering the neutral birth- 
death model. 

Solution of the model - We study the master equation of 
the stochastic process describing the non-neutral dynamics 
taking place in the model ecosystem. If we call N n (f, t) the 
number of species of fitness / populated by n individuals, 
its master equation is given by 



dN n (f,t) 
dt 



yN n _ 1 (f,t)-b n N n (f,t) 



d n+1 N n+1 (f, t) - d n N n (f, t) + Np{f)S nil , (8) 
where the birth rates b n and death rates d n are given by 

K(f) = £*(/)-(£) 2 C{f) 

dn(f) = ^D(f)-Qfc(f) (9) 

with B(f) = ae 0S yZi + 7, D(f) = 1 + a + ^e~ 0S / Z 2 and 
C(f) = ae 13 ! / Zi + je^Pf /Z 2 and where we have defined 
Zi , Z 2 as 

Zi(t) = 1 [ dfY,^ h nN n {f,t) (10) 
Z 2 (t) = j- f dfJ2^ f *nN n (f,t)- 

J n 

In order to solve the master equation (8) we assume that 
Z\{t) and Z 2 (t) in the asymptotic limit t — > 00 converge 
to a constant value 



Z 1/2 (t)^Z 1/2 + 0(N~ 1 / 2 ) 



(11) 



In this assumption, the steady state of the master equation 
(8) can be solved exactly as in neutral birth-death models 
[4,5] giving the distribution of the abundance of the species 
with given fitness. The stationary solution for N n (f) is 
given by 



Nn(f) = N 



P(.f) fr bi-iif) 



(12) 



In the limit in which n <C N this distribution becomes a 
convolution of different Fisher log series distribution [5] 



N n (f) = N 



p(f) W)} 

D(f) 



n-l 



(13) 



with the diversity number 9 — 9(f) dependent on the fit- 
ness of the species / according to the relation 



e(f) = B(f)/D(f) 



7 



1 e~Pf/Z 2 + (a + 1) ' 



(14) 



From now on we focus on this model in the case of 
bounded fitness distributions with / £ [/ m , /m] for which 
we might expect that, al least in a certain region of the 
phase space we might have a stationary state of this evo- 
lutionary dynamics. Without loss of generality we take 
from here on / <E [0, 1]. In this case the diversity parame- 
ter 9(f) reaches a maximum at the maximal fitness / = 1. 
The average number of individuals with a given fitness / 
is 



<n(/)> = N 



P(f) 



1 + a - 7 + 7e-W/Z 2 - aePf/Zi 



(15) 



This distribution admits two relevant limits a = and 
7 = 0. To show this, it is useful to put / = 1 — e and to 
define a distribution g(e) — p(l—e). In absence of negative 
selection, i.e. 7 = 0, the average number of individuals of 
fitness / = 1 — £, (n(f)) can be expressed in terms of a 
Bose-Einstein distribution of the occupation of the energy 
level e. 

Conversely when reproduction selection are absent, i.e. 
a = 0, the average number of species (n(f)) with fitness 
/ = 1 — e follows a Fermi-Dirac occupation of the energy 
levels e, provided that the negative selection processes oc- 
cur at a smaller rate than speciation, i.e. when 7 < 1. 
This show clearly that the distribution of species in this 
evolutionary dynamics change significantly if only the re- 
production selection process occurs or if only to the neg- 
ative selection process occurs. Moreover taking the limit 
[3 — > 00 for a = 0, 7 < 1, we recover for (n(f)) the step 
function characteristic of the Bak-Sncppen model [14] and 
we can map the dynamics to the dynamics of a mean-field 
Bak-Sncppen model [15] showing punctuated equilibrium. 

Eq. (13) and Eq. (15) are well defined once Z\ and Z 2 
are found self-consistcntly. The self-consistent equations 
read 

Zl = J 1 + a - 7 + ieP*/Z 2 - ae^/Z, 



= Jo dEg ^l 



+ a — 7 + 7e^ e I Z 2 — ae @ e / Z\ 



(16) 



These self-consistent equations have to be solved for 
each given parameter value (a, 7,/?) and each distribu- 
tion of the fitness p(f) = g(e = 1 — /). Eqs. (16) might 
lack a solution when the maximal fitness value is reached 
in the system at infinite time, i.e. when lim/^/ M p( f ) = 
\mi e ^o g(s) = 0. In this case, a inverse critical tempera- 
ture can be evaluated [21] such that for (i < [3 C the Eqs. 
(16) have a solution while for (3 > (3 C they do not. This 
phase transition can be mapped in the case 7 = to 
a Bose-Einstein condensation phase transition in a Bose 
gas. As in the Bose-Einstein condensation the lowest en- 
ergy state gets a finite occupation, in this model when 
(3 < (3 C we observe the emergence of dominant species in 
the ecosystem having a finite fraction of the total popula- 
tion. This dominant species can change with time and the 
ecological system is non-stationary. We observe here that 
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Fig. 1: Fraction of the individuals of the population with fitness 
/ = 1 — s for a uniform distribution p(f) with / £ [0, 1]. The 
data are shown for f3 = 20 and are averaged over 100 runs of 
the non-neutral birth-death dynamics involving a population 
oi N = 1000 individuals. The data are collected at the time 
T — 5 x 10 4 of the stochastic dynamics. 



the time needed for an increase in the average fitness of 
the ecosystem in the condensed state increase with time 
in a multiplicative way. In fact we can say that this evo- 
lutionary model show in this limit relevant aging effects. 
This aging phenomena might be quite unrealistic for real 
ecologies. This feature of the model can be overcome by 
including genetic mutations that can occur on very long 
time scales. Genetic mutations are therefore negligible if 
we consider an ecosystem which reach quickly the steady 
states, but they can be have relevant effects on the time 
scales of the aging dynamics of the ecosystem in the con- 
densed phase. 

Numerical characterization of the stationary and non- 
stationary state of the ecosystem - We have simulated 
the model and studied the population distribution and 
collected numerical evidence for the condensation phase 
transition to a non-stationary state. In fig. 1 we report 
the total population (n(f)) as a function of the fitness 
/ = 1 — e for a uniform fitness distribution p(f) = 1 
with / e (0,1). The figure shows that at low tempera- 
ture, in absence of the reproduction process, high fitness 
species are not rewarded by the non-neutral birth-death 
dynamics. On the other hand, in absence of the negative 
selection process, low fitness species are not punished by 
the non neutral birth-death dynamics. In order to study 
the non-stationary condensed state, we consider a fitness 
distributions of the type p(f) = (k + 1)(1 — f) K . 

In fig. 2 we report the time dependence of the fraction 
of individuals in the most populated species above and be- 
low the condensation transition. Above the condensation 
phase transition only an infinitesimal fraction of the total 
population belongs to the most abundant species, while 
below the condensation transition the dominant species 
is populated by a finite fraction of all individuals. The 
dominant species changes over long time scales as it is 
eventually overcome by other fitter species arriving later 



Fig. 2: Time dependence of the fraction of individuals in the 
dominating species below (/3 = 10) and above (/3 — 0.1) the 
condensation phase transition. Below the phase transition the 
non-stationary of the process is manifest with the most domi- 
nating species changing over long time scales. The data repre- 
sent a single run of the evolutionary dynamics with parameters 
N = 1000, a = 10 and 7 = 0.1 and k = 0.1. The critical value 
of the inverse temperature for these parameters is given by 
/3c = 2.4. 



in the ecosystem. This process is the process that gener- 
alize the principle of competitive exclusion to open ecolo- 
gies. In close ecolgies the competitive exclusion principle 
in fact states that at long time scales, the fittest species 
of the ecosystem is fixated with higher probability. Here 
we prove that, in the presence of migration, the ecology 
might become dominated by a species with high fitness. 

In order to characterize the phase transition as a func- 
tion of (3 we studied the order parameter n max /N indicat- 
ing the fraction of individuals in the dominating species. 
Going across the condensation phase transition toward the 
non-stationary condensed state, the diversity of the pop- 
ulation rapidly decreases as it is shown in Fig. 3 where 
we plot the total number of species S/N as a function of 
(3. Moreover the inhomogeneity in the systems increases 
strongly below the condensation phase transition. In or- 
der to measure this inhomogeneity we plot the parameter 8 
defined as 5 = — \'^g 2 ^ with Y2 indicating the participation 

ratio of the species population, i.e. Y% = J2i=i (lv) • 

The participation Y2 varies from Y2 = l/S when all the 
species are equally populated to Y2 = 1 when one species 
dominates the population. Therefore the parameter S G 
(0, 1) is equal to S = 1 for equally populated species and is 
5 = when only one species is occupied. In fig. 3 we show 
how n max /N, S/N, and S varies across the condensation 
transition as a function of (3. 

Conclusions - In this paper we have proposed a new, 
simple, non-neutral stochastic model for ecosystems in 
which species compete for finite resources. In the model 
we distinguish different limiting cases: 

• Case (3 = - We recover the Neutral Model of the Uni- 
fied Theory of Biodiversity with species abundance 
distribution that can be approximated by a log Fisher 
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Fig. 3: The normalized number of species S/N, the fraction 
of individuals in the dominating species n max /N and the pa- 
rameter 5 which measures the diversity in the population as a 
function of /3 across the phase transition for a = 1.5, 7 = 1.5 
and k = 1.0, N = 1000. The solid line indicates the value of the 
critical inverse temperature. 



distribution with diversity parameter given by Eq. 
(14), therefore 6= (a + j)/(l + a + 7). 

• Case (3 > j3 c - We have a stationary ecology where the 
species abundance can be approximated by a convo- 
lution of log Fisher distributions with diversity pa- 
rameter dependent on the fitness 9 = 0(f) given by 
Eq. (14) and constants Z\ and Zi satisfying Eq.(16). 

• Case f3 < (3 C - We have a non stationary ecology with 
a dominating species and strongly reduced diversity. 
The phase transition at (3 = j3 c extends the princi- 
ple of competitive exclusion to open ecosystems with 
high migration rate and shed light on the instabilities 
found in presence of high fitness invasive species. 

• Case a = and j3 — * 00 - We find a dynamics of 
extinctions that can be mapped to a mean-field Bak- 
Sncppen model and it can be shown [21] that this 
model shows punctuated equilibrium. 

We believe that the proposed model can provide a gen- 
eral framework for the study of non-neutral models of 
biodiversity presenting in a unified framework different 
known limits of stochastic dynamics of ecosystems. Work 
on progress is investigating extensions of this model for the 
description of ecosystems with non-trivial foodweb prey- 
predator or mutualistic interactions. 
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